#ifndef EQUATIONSOLVER
#define EQUATIONSOLVER

#include "Function.hpp"
#include <cmath>
#include <iostream>
using namespace std;
class EquationSolver{
	protected:
		const Function & F;
	public:
		EquationSolver(const Function& F): F(F) {}
		virtual double solve()=0;

};
class Bisection_Method : public EquationSolver{
	private:
		double a,b;
		double eps,delta;
		int Maxiter;
	public:
		Bisection_Method(const Function &F,double a,double b,
				double eps=1e-7,double delta=1e-6,int Maxiter=50):
				EquationSolver(F),a(a),b(b),eps(eps),delta(delta),Maxiter(Maxiter){}
		int sgn(double x){
			return (x>0)?1:(x<0)?-1:0;
		}
		virtual double solve(){
			double u=F(a);
			double v=F(b);
			double c,w;
			for(int i=0;i<Maxiter;++i){
				double h=b-a;
				c=a+h/2;
				if(abs(h)<delta){
					break;
				}
				w=F(c);
				if(abs(w)<eps){
					break;
				}else if(sgn(w)!=sgn(u)){
					b=c;
					v=w;
				}else{
					a=c;
					u=w;
				}
			}
			return c;
		}
};
class Newton_Method:public EquationSolver{
	private:
		double x0;
		double eps;
		int Maxiter;
	public:
		Newton_Method(const Function &F,double x0,
				double eps=1e-7,int Maxiter=50):
			EquationSolver(F),x0(x0),eps(eps),Maxiter(Maxiter){}
		virtual double solve(){
			double x=x0;
			for(int i=0;i<Maxiter;++i){
				double u=F(x);
				if(abs(u)<eps){break;}
				x=x-u/F.derivative(x);
			}
			return x;
		}
};
class Secant_Method: public EquationSolver{
	private:
		double eps,delta;
		int Maxiter;
		double x_0,x_1;
	public:
		Secant_Method(const Function &F,double x_0,double x_1,
				double eps=1e-7,double delta=1e-6,int Maxiter=50):
                        EquationSolver(F),x_0(x_0),x_1(x_1),eps(eps),delta(delta),Maxiter(Maxiter){}
		virtual double solve(){
			double x_n=x_1;
			double x_n_1=x_0;
			double u=F(x_n);
			double v=F(x_n_1);
			for(int i=2;i<=Maxiter;++i){
				double s=(x_n-x_n_1)/(u-v);
				x_n_1=x_n;
				v=u;
				x_n=x_n-u*s;
				if(abs(x_n-x_n_1)<delta){
					break;
				}
				u=F(x_n);
				if(abs(u)<eps){break;}
			}
			return x_n;
		}
};
#endif
